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ABSTRACT 

Radiation feedback from young star clusters embedded in giant molecular clouds (GMCs) is believed to 
be important to the control of star formation. For the most massive and dense clouds, including those 
in which super star clusters (SSCs) are born, pressure from reprocessed radiation exerted on dust grains 
may disperse a significant portion of the cloud mass back into the interstellar medium (ISM). Using 
our radiation hydrodynamics (RHD) code, Hyperion, we conduct a series of numerical simulations 
to test this idea. Our models follow the evolution of self-gravitating, strongly turbulent clouds in 
which collapsing regions are replaced by radiating sink particles representing stellar clusters. We 
evaluate the dependence of the star formation efficiency (SFE) on the size and mass of the cloud 
and k, the opacity of the gas to infrared (IR) radiation. We find that the single most important 
parameter determining the evolutionary outcome is k, with n ^ 15 cm 2 g -1 needed to disrupt clouds. 

For k = 20 — 40 cm 2 g -1 , the resulting SFE = 50 — 70% is similar to empirical estimates for some 
SSC-forming clouds. The opacities required for GMC disruption likely apply only in dust-enriched 
environments. We find that the subgrid model approach of boosting the direct radiation force L/c by 
a “trapping factor” equal to a cloud’s mean IR optical depth can overestimate the true radiation force 
by factors of ~ 4 — 5. We conclude that feedback from reprocessed IR radiation alone is unlikely to 
significantly reduce star formation within GMCs unless their dust abundances or cluster light-to-mass 
ratios are enhanced. 

Keywords: hydrodynamics - methods: numerical - radiation: dynamics - radiative transfer ISM: 
clouds - stars: formation - galaxies:star clusters 


1. INTRODUCTION 

Giant molecular clouds (GMCs), the sites of star for¬ 
mation, form out of the diffuse interstellar medium (ISM) 
due to some combination of self-gravity, the gravity of the 
stellar disk and bulge (which compresses the ISM verti¬ 
cally everywhere and horizontally in spiral arms), and 
large-scale gas motions associated with turbulence, su- 
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nent structures, and based on age-dating of associated 
star clusters (in the Milky Way and other galaxies) are 
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2014). It is widely be- 


for the demise of GMCs, but exactly how this works is 
still poorly understood. Among the most basic uncer¬ 
tainties is which among the possible feedback effects pre¬ 
dominate for different regimes of GMC properties and 
surrounding environment. Currently, the most-discussed 
can didate effects (see, e.g., the review of Krumholz et al. 
20141 include (1) (magneto)hydrodynamic" forces from 
overpressured regions produced by ionizing radiation, 
shocked stellar winds, or supernova blasts; and (2) radia¬ 
tion forces from the primary absorption of stellar optical 
and ultraviolet (UV), and from the secondary absorption 
of reprocessed infrared (IR). 

Both theoretical and observational motivations have 
led to an increased interest in the effects of radiation 
forces. Although supernovae inject an order of magni¬ 


tude more momentum per stellar mas s to their surro und¬ 
ings t han other forms of feedback (Ostriker fe S hctty 
201 l|j |Kim fe Ostriker |2015| |lttrig &; ilennebelle||2()15| 


Walch & Naab 2014 Martizzi et al. 2015[ (Teen' et ah] 

2015|), there is a significant time delay between the ad- 


vent of star formation in a GMC and the explosion of the 
first supernova (and further delay before the last massive 
stars die). If other agents are able to destroy the GMC 
in this interval, much of the momentum and energy from 
supernovae (SNe) may be delivered to the diffuse ISM 
rather than the progenitor’s birth cloud. 

In GMCs with low escape speeds hosting clusters con¬ 
taining massive stars, the combination of photoevapo¬ 
ration and the pressure force from the expanding H II 
region can unbind much of a cloud’s initial mass (e.g., 


therein). However, analytic spherica 

models (Krumholz 

& Matzner 2009 

Murray et al. 
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clusters within clouds suggest t 
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le effects from radi- 


ation pressure will exceed that from ionized gas pressure 
at high values of the total cloud mass and surface density. 
Furthermore, the relative importance of reprocessed IR 
compared to direct optical/UV is expected to increase 
as the cloud surface density increases. Observational ev¬ 
idence suggests that radiation pressure exceeds ionized 
gas pressure close to the centers of H II regions around 
massiv e clusters (consistent wit h the theory of |Draine 


2011a), and for youn ger systems (Lopez et al 12011 2014 


Pellegrini et al.|2007 1. Although historically the gas pres- 
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sure from shocked stellar winds was expected to drive dy¬ 
namics around massive clusters at early times, this has 
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been called into question for some systems due to the lack 
of X-ray emission, with the suggestion that th e primary 
wind and/or hot shocked gas largely escapes ( Townslcy 


et al.||2003[ |Harper-Clark &; Murr ay 2009; | Rogers & Pit- 

tard||2013|j. 


Embedded super star clusters (SSCs) in dense GMCs 
represent the systems for which radiation pressure is ex¬ 
pected to play the greatest role. Considering that they 
are still deeply buried in GMCs, the primary signature 
of these SSCs is thermal radio emission, indicating clus¬ 
ters of masses ~ 10 4 — 10 6 Mq powering H II regions 


only a few pc in size 

(Turner et al. 1998, 2000 Kob- 

ulnicky & Johnson 

1999; Johnson e 

al. 2001; Johnson 

& Kobulnicky 2003 

J 

ohnson et al. 

20151 Reines et al. 

2008 Tsai et al. 2009 

Kepley et al. 
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2014). The molec- 


than typical Milky Way GMCs. For example, “Cloud 
D” in NGC 5253 has a mass ~ 2 x 10 6 M 0 and diame¬ 
ter ~ 40 pc based on the SMA observations of |Turner| 


et al. (2015), giving a mean density of hydrogen nu¬ 
clei njj = 1800 cm -3 and mean surface density X = 
1600 Mq pc -2 = 0.33 g cm -3 . The “pre-SSC” cloud 


in the Ant ennae system studied with ALMA by |John- 
son et al. (2015) has an estimated mass ~ 0.3 — 1.5 x 
10 Y Mq and diameter < 40 pc, yielding nu ~ 0.26 — 
1.3 x 10 4 cm" 3 and X > 0.2 - 1.2 x 10 4 M 0 pc" 2 . In 
the central starburst region of NGC 253, ALMA ob ser¬ 
vation s of high critical density tracers by [Leroy et al. 


(2015) identified clouds with radii 10 — 50 pc and masses 
0.2 — 6 x 10 7 M 0 , implying typical nn ~ 2000 cm -3 and 
X ~ 6000 Mq pc" 2 . 

The deeply embedded nature of SSCs makes it difficult 
to constrain the relative importance of different feedback 
mechanisms empirically. However, the SFEs appear to 
be at least an order of m agnitud e hig her t han in typical 


Milky Way GMCs (e.g., Tsai et al. 2009). The young, 


isolated Cloud D in NGC 5253 has an estimated star 
formation efficiency (SFE) ~ 0.6, well above the range 
0.01 — 0.2 estimated fo r the most lum inous Milky Way 


star-forming complexes (Murray 2011). It is clearly of 
interest to develop theoretical models that explore the 
effects of feedback in controlling star formation and set¬ 
ting SFEs in GMCs comparable to SSC hosts. This may 
also inform our understanding of globular cluster forma¬ 
tion, which presumably occurred in similarly dense and 
massive clouds. 

Because forces from reprocessed IR radiation become 
increasingly important as gas surface densities increase^ 
the pressure from trapped IR has been argued to domi- 


systems ( 

Thompson et al. 2005| 

Murray et al. 

2010) 

However, analytic models typicalJ 

y require major sim 


plifications to be tractable (including having all the gas 
mass collected in a thin, uniform-density spherical shell, 
and having a single stellar cluster that is centrally located 
with a luminosity that is fixed in time), and it is unclear 
the extent to which these may affect the conclusions. It 
is therefore useful to employ time-dependent numerical 
simulations to investigate systems with more structural 

1 For the spherical case, the total force imparted to the gas 
is (1 + tir)L,/c where L* is the stellar luminosity and qq is the 
center-to-edge IR optical depth, proportional to the cloud’s surface 
density. 


and temporal complexity. 

In this paper, we shall consider numerical models of 
massive, compact, turbulent GMCs that fragment grav¬ 
itationally to form massive star clusters. Our model 
clouds have surface density 1600 — 6200 M 0 pc -2 , com¬ 
parable to the GMCs in starburst nuclei within which 
SSCs are born. Optical and UV radiation from young, 
hot stars dominate the luminosity of massive clusters, 
but these primary photons are likely to be absorbed by 
dust very near their source. This warm dust then ra¬ 
diates isotropically into the lower-frequency IR band, 
which has a much smaller absorption cross-section but 
can still be absorbed and re-emitted multiple times in 
high-column GMCs. 

In our simulations, we apply radiation hydrodynam¬ 
ics (RHD) methods to focus on the effects of long- 
wavelength radiation. We consider the regime of clouds 
at relatively high mass and density, such that the optical 
depth to IR exceeds unity and simple spherical models 
would predict that the radiation force from reprocessed 
IR exceeds that of the direct optical/UV. We follow the 
co-evolution of the gas and radiation in the system, under 
the assumption that radiation forces applied to dust are 
transferred to the gas (i.e., assuming perfect gas-dust col- 
lisional coupling), and that all radiation that is absorbed 
is locally re-emitted. Over time, the fraction of a cloud’s 
mass converted to stars increases as more and more gas 
collapses. However, as the stellar luminosity grows, the 
radiation field exerts increasing outward forces on the 
gas. Under certain conditions, we show that a significant 
portion of the cloud’s initial gas mass is expelled from 
the system. 

Our models are highly idealized, in that we con¬ 
sider exclusively the effects of long-wavelength radiation, 
with spatially-uniform opacity. As our numerical code 
presently allows only a single opacity, here we do not con¬ 
sider the direct effects of either ionizing or non-ionizing 
UV; effectively, we treat radiation as being degraded to 
IR close to each stellar cluster source. In reality, several 
feedback effects operate simultaneously, and complete 
models must eventually include (at least) non-ionizing as 
well as ionizing UV, and stellar winds. The present set of 
simulations provides a baseline for more comprehensive 
studies, similar to the baseline provided by simulations 
that focus exclusively on the effects of ionizing radiation 
and cloud disruption from H II region e xpansion (e.g., 
Dale et ah|2012 2013 Walch et al.||2012 ). 


We note that all of our clouds have escape speeds ex- 
ceeding 20 km s - 1 . In this regime, the simulations of 


Dale et al. (2012 2013) suggest that photoevaporation 
combined with the pressure from ionized gas would be 
able to unbind no more than 1% of the cloud’s mass 
prior to the advent of supernovae. Additionally, al¬ 
though realistic RHD models focused on non-ionizing 
UV radiation have not yet been completed for this ex¬ 
treme regime (very high X and large v esc ), analytic mod¬ 
els have been developed that account for the effects of 
turbulence- driven internal structure in the radiation/gas 
interaction (Thompson & Krumholz| 2014; Raskutti, Os¬ 
triker, & Skinner 2015, in preparation); the Raskutti et 
al model has also been verified via numerical RHD simu¬ 
lations in the lower-X regime. These models predict that 
it would be difficult for the direct UV radiation to ex- 
pell substantial material from large, dense GMCs, since 


















































































RHD Simulations of Turbulent, Star-Forming Clouds 


3 


the already-high mean surface density is increased in the 
filamentary structures that comprise most of the cloud’s 
mass. Thus, as prior studies suggest that neither ion¬ 
izing or non-ionizing UV will be effective in destroying 
GMCs for the cloud regime studied here, we are moti¬ 
vated to turn the focus to the effects of reprocessed IR. 
To our knowledge, this work represents the first direct 
RHD study of the dynamical effects of reprocessed radi¬ 
ation in turbulent, self-gravitating, star-forming clouds. 

The plan of this paper is as follows. We begin with 
a summary of our numerical methods, model specifica¬ 
tion, and model parameters^together with tests to verify 
code performance (Section [ 2 ]). We th en consider evolu¬ 
tion of a fiducial model (Section 3.1), which also serves 


to illustrate the radiation structure and differences in the 
radiation/matter inte racti on compared to s imple spher¬ 
ical systems (Section 3.2). In Section 3.3 we analyze 
the effects of varying the opacity, demonstrating that 
quite high opacity (re > 10 cm 2 g _1 ) is required for sub¬ 
stantial mass to be ejected by radiation forces; we also 
quantify the gas-radiation anticorrelation and the level 
of radiation trapping. Our set of simulations allows for 
varying clo ud m ass and radius as well as a range of opac¬ 
ity; Section [3 . 4| compares simulation outcomes (including 
net SFE and net momentum ejected) for different model 
parameters. Section[4]summarizes and discusses our con¬ 
clusions. 

2. NUMERICAL METHODS & MODEL DESCRIPTIONS 


2.1. Numerical Methods 
We evolve the equatio ns of RHD usin g 


duno v co de Hyperion (Skinner & Ostriker 2013 
after S013 1. Hype rion is an extension of the Athena code 


our Go- 
here- 


(Stone et ah 2008) for computational hydrodynamics and 
magnetohydr odynamics (MHD); we em ploy the van Leer 
algorithm of Stone & Gardiner (2009) to integrate the 
gas equations, utilizing the HLLC Riemann solver and a 
piecewise-linear spatial reconstruction scheme. For sim¬ 
plicity in this first study, we neglect magnetic fields and 
adopt an isothermal equation of state (EOS); as the 
sound speed is small compared to other speeds in the 
problem, the dynamics are insensitive to its exact value. 
Gravity of the gas as well as the “star particles” (repre¬ 
senting collapsed regions that have formed star clusters) 
is computed via Fourier methods, as described below. 

In Hyperion, the radiation energy density and flux, i.e., 
the zeroth and first moments of the radiation intensity, 
are advanced in time using a piecewise-linear spatial re¬ 
construction and an HLL-type Riemann solver. For the 
radiation energy equation, we adopt the limit of radia¬ 
tive equilibrium, such that all absorbed radiation is re¬ 
emitted locally. For long wavelength radiation, which we 
consider, scattering is small compared to true absorp¬ 
tion and is neglected in our treatment. As velocities are 
small compared to the speed of light and optical depths 
are moderate, we adopt the static diffusion limit in which 
terms of 0{v/c) and Oirv/c) are neglected. Finally, we 
employ the reduced spe ed of light approximation (RSLA) 
(Gnedin & Abel 2001), which allows us to solve the ra- 
diat ion subsystem explicitly rather than implicitly (see 
also Rosdahl et al. 2013). The only sources of radiation 
are the star particles. 

The system of equations to be solved for the gas and 


radiation is given by 

d t p + V • (pv) =0, 


(la) 


<9 t (pv) + V • (pvv + PI) = -pV$ + pre —, (lb) 

c 


IdtS + V- (-) = -, 

c \c J c 

1 /F\ F 

- d t — + V • P = —pre —, 
c V c / c 


(lc) 

(l d) 


where p, v, and P are the gas density, velocity, and pres¬ 
sure, respectively, and $ is the gravitational potential. In 
Equations ([T]), £, F, and P are the frequency-integrated 
radiation energy density, flux vector, and pressure tensor, 
respectively, which are measured in the inertial frame, 
and c is the reduced speed of light. The absorption opac¬ 
ity, re, is taken to be a spatial constant, although we vary 
this parameter in different models. Finally, the term j* 
in Equation flc| ) represents emission of radiation from 
star particles. 


The Hyperion code uses the Mi closure relat ion (Lever- 
more & Pomraning|1981 Gonzalez et al.|2007 ) to express 
P as a function of £ and F in Equations (lj. The radi¬ 
ation subsystem in Equations (lc) and (Id ris operator- 
split from the gas subsystem in Liquations | la[) and © 
and is explicitly evolved by subcycling on a time step 
determined by the radiation propagation speed, c. Us¬ 
ing the RSLA, we choos e c satisfying the RSLA static 
diffusion criterion (S013), 

c> 

^max max{ 1, T max }, (2) 

where u max is the maximum velocity of the gas and r max 
is the maximum optical depth in the simulation. This al¬ 
lows us to reduce the radiation-to-gas time step ratio to 
a reasonable level, but has a negligible effect on the sys¬ 
tem’s dynamics as the gas evolves much more slowly than 
the radiation field. Our simulations typically use a value 
of c such that c/u max ~ 10 — 100, which is substantially 
less than realistic values of c/i; max ~ 10 4 — 10 5 , making 
explicit time integration computationally feasible. 

In Hyperion, star particles are treated us ing the algo¬ 
rithm developed by Gong & Ostriker (20131, with added 
functionality to provide for the luminosity of the parti¬ 
cles)^] The Mi closure cannot resolve the behavior of a 
streaming radiation field too near to a true point source, 
since the the radiation flux from such a source varies 
rapidly in angle. Therefore, radiation sources in our al¬ 
gorithm must be resolved over some minimum number 
of grid zones. We have found it convenient to add radi¬ 
ation energy density to the grid using a Gaussian source 
function given by 


j*( x ) = 


( 2 w 2 ) 3 / 2 


exp - 


2a 2 


(3) 


where L* is the star particle’s luminosity, x* is the star 
particle’s position, and cr* = R*/V 2 In 2 is set such that 
the half-width at half-maximum (HWHM) of the distri¬ 
bution is equal to the star particle’s effective size, i?*. 
Note that / Q r j*(r') 47rr' 2 dr' —> L* for r/R* 1. In 

2 At the resolution of our simulations, each “star particle” actu¬ 
ally represents a fully-sampled star cluster. 
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practice, we have found that sources with R*/Ax > 8 
are sufficiently well-resolved that angular variations in 
the rad iatio n flux at radii r i?* are negligible (see 
Section |2(2|). 

To compute gravitational forces including contribu¬ 
tions from both gas and star particles, we first use the 
particle mesh (PM) method to assign the star particle 
masses to a discrete grid via th e triangular-shaped clou d 
(TSC) method, as described in Gong fe Ostrikeij (2013 ). 
We then apply the “zero-padding” method of |Hockney| 
& Eastwood (1988) to obtain the potential, <f>(x), of an 
isolated source distribution subject to open (vacuum) 
boundary conditions via fast Fourier transforms (FFTs). 
This potential is given by the solution of Poisson’s equa¬ 
tion, 

V 2 <fi = 47rGp(x), (4) 

for a given density field, p(x); p includes both contribu¬ 
tions from the gas and the star particles. Solutions of 
Equation Q can be expressed as the convolution 


4>(x) = G /S(x,x>(x')d 3 x', 


(5) 


where C/(x,x') = C?(|x — x'|) = — |x — x'| 1 is a Green 
function solution of the equation V 2 t/ = 47rd 3 (x — x')- 
Using the Fourier convolution theorem, Equation (|5j) 
may be re-expressed as the convolution of the Fourier 
transforms of the density distribution and Green func¬ 
tion; details of this are given in the Appendix. 


2.2. Model Description 

Each model cloud is initiated as a uniform-density 
sphere of radius i? c ioud and mass M c i ou( j. The cloud is 
centered in a computational box of side length Lbox = 
4i?doud with background density set to 1% of the ini¬ 
tial cloud densityrj A highly supersonic turbulent veloc¬ 
ity field is appliecfinitially, which rapidly creates density 
structure within the cloud. Over time, this initial veloc¬ 
ity field decays, although collapse and radiation forces 
drive further turbulent motions. 

We initialize the turbulent velocity field using Gaussian 
random perturbations with power spectrum |<5v| ex fc~ 4 , 
for k/dk g [2 64], where dk = 27r/Lbox, as described 


Stone et al. (1998). The velocity perturbations are 


normalized such that a v j rj i n j t = 2, corresponding to a 
just-bound state, where 


= 2 


Akin 


E, 


( 6 ) 


grav 


is the virial parameter, and where -Ekin = (l/2)M c i ou d^' 2 
and E grav = (3/5)GM c i ou d 2 /-R c ioud are the initial ki¬ 
netic and gravitational energies of the gas, respectively; 
a = [a v i r? i n it3GM c i ou d/(5i?cioud)] is the initial turbulent 
velocity dispersion^] The velocities of the perturbations 
are recentered such that they add no net momentum to 
the computational domain, i.e., such that f pSv dV = 0. 


3 Realistic clouds may have lower density contrast relative to 
their surroundings, which could lead to additional late-time accre¬ 
tion. However, for this first study we consider isolated clouds for 
simplicity. 

4 Note that for a uniformly dense sphere, such as our initial 
configuration, a v i r is equivalent to the often-us ed definition of the 
virial parameter from | Bertoldi &; Mc Kee (1992 , see Equation 2.8a). 


For all simulations, the isothermal sound speed is set 
to c s = 2 km s -1 . We adopt this value, somewhat larger 
than true thermal sound speeds in cold clouds, in part 
to limit extreme shock compression, as would occur if 
we had included magnetic fields. We have found that 
results are insensitive to the exact choice of the sound 
speed, provided the Mach number is large. The initial 
turbulent Mach number is ~ 11 for our fiducial parame¬ 
ters (see below), and more generally the turbulent flows 
in all of our models are highly supersonic. Tests show 
that there is little variation in simulation outcomes for 
different realizations o f the initial random perturbed ve¬ 
locity field (see Section 3.1), so we use a single realization 


for most set s of model parameters. 

Following Gong & Ostriker (2013), the density thresh¬ 
old for star particle creation is set using the Larson- 
Penston criterion, p t hr = Ppp(t = Ax/2 ), where 


Php{r) = 


8.86c 2 

47rGr 2 


(7) 


describes the density profile of the asymptotic state pro- 
duced by gravitational collap se of an isothermal sphere 
(Larson 1969 Penston 19691; this profile is approached 
tor collapse initiated from a wide variety of initi al condi¬ 


tions, including supe rsonic turbulent f lows (see Gong fe 
Qstrik er||2009l |2011|). As discussed in |Gong & Ostriker| 
j 2013|, tne Larson-Penston density criterion is a factor of 


14 times larger than the Truelove criterion (Truelove 
et al.|1997 ) used in man y other star particle cr eation lm- 


et al.|i99f I used m man y otner star particle creation im¬ 
plementations (s ee, e.g.,|Krumholz et al.|2004 Federrath 
& Klessen||2012|), but the star particles produced in a 
Hi 


particles produced 
turbulent How follow essentially the same histories. 

We assume the star clusters that form in the cloud 
(represented as star particles) have a bolometric luminos¬ 
ity per unit mass of T = 17 00 er g s~ 4 g ~ 4 , as o btained 


from a Starburst99 model ( Leitherer et al.|fl999 ) for the 
total luminosity of a young cluster that fully samples a 
Kroupa initial mass function (IMF), averaged over a life¬ 
time of 5 Myr(^] Since the total luminosity of a cluster 
would not change substantially over the cloud’s lifetime, 
for simplicity we set L* = 'PM*, independent of the star 
particle’s age. 

We adopt a constant eHective physical size for each 
star particle, which we take to be f?* = 1 pc. Sources 
of this size are generally consistent with obs ervations of 


young, embedded super-star clusters (e.g., Johnson & 
Kobulnicky 2003) as well as older clusters ( |Murray|2009p . 


INote that here, we do not attempt to model the H II 
region that those star clusters would create. For high lu¬ 
minosity sources, UV radiation pressure on dust signifi¬ 
cantly alters the conditions within H II regions compared 


to the classical Stromgren solutio n (Binette et al. 
Dopita et al.||2002 Draine 2011a), with radiation 


1997 


orces 


strongly compressing the photoionized gas. In principle, 
high pressure from a shocked stellar wind could further 


5 For random sampling of a fixed IMF, IP saturates for clusters 
of mass M c \ > 10 4 Mq. In our simulations, we find that star 
particles, which represent star clusters on the scales we consider, 
may occasionally form with masses of order a few times less than 
10 4 Mq. However, these somewhat under-sampled clusters typi¬ 
cally comprise a trivial fraction of the total mass in stars, which 
is dominated by fully-sampled clusters. Furthermore, these lower- 
mass star particles tend to accrete sufficient gas over a short enough 
time scale to justify the use of a constant specific luminosity. 
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flux-, and Rosseland-mean opacities, which are defined, 
respectively, by 


Figure 1. Time evolution of u max Tgo% for three different values 
of the reduced speed of light, c, where u max = a c s approximates 
the maximum hydrodynamic signal speed and Tg 0 % is the optical 
depth at the 90th percentile over all projections of the grid in all 
directions. All runs were performed with N = 128 3 zones. This 
demonstrates that the time evolution of the optical depth for the 
bulk of the gas is insensitive to the actual value of the reduced 
speed of light used, provided the RSLA static diffusion criterion 
of Equation pi is satisfied. The value of c = 4900 km s —1 is the 
reduced speea of light in our fiducial run. 


compress the photoionized gas, but observed ionization 
paramet er measures suggest th at this does not occur in 
practice (Yeh & Matzner 2012). Given the uncertainties 


about the effective size of the emission region in which 
the cluster’s radiation is reprocessed into IR, we opt to 
keep i?* fixed rather than having it depend on the mass 
and age of each star particle. Tests indicate that varying 
the physical size i?* of each star particle radiation source 
by a factor of two in either direction at fixed grid reso¬ 
lution only affec ts t he final total mass in stars by about 
5% (see Section |3.1 ). 

From Equation ([2]), the required value of c depends on 
the maximum optical depth and velocity in the model. 
If we define rg 0 % at a given time as the optical depth at 
the 90th percentile over all projections of the grid in all 
directions, our tests show that the maximum of Tg 0 % over 
a simulation run is comparable to r, the initial IR optical 
depth across the GMC, within a factor of a few, and 
that this value is well-converged if we use c = 5u max r. 
Figure [l] shows the time evolution of u max Tgo% for three 
different values of c, where u max = a + c s approximates 
the maximum hydrodynamic signal speed. Each of these 
runs was performed at a reduced resolution of N = 128 3 
zones for parameters otherwise the same as the fiducial 
model (K20), and the choices of 2500, 4900, and 9900 
km s -1 for c correspond to 5, 10, and 20 times (<j + c s )t, 
respectively. The data clearly indicate that the optical 
depth that is “seen” by the majority of the radiation is 
insensitive to the specific value of c, provided the RSLA 
static diffusion criterion of Equation ([2]) is satisfied. In 
all of our simulations, we therefore choose c such that c > 
5(cr + c s )r; in most simulations we have c = 10(cr + c s )t. 

In our current investigation, which focuses on effects of 
reprocessed radiation that is absorbed and re-emitted by 
dust, the absorption opacity of the medium, k , is taken 
to be a constant in space and time. More realistically, 
the opacity law would depend on the frequency of the 
radiation as well as on local dust abundance and prop¬ 
erties, such that the frequency-averaged opacity would 
depend on these properties and also the particular radi¬ 
ation regime. For example, consider the energy-, Planck-, 


_ / du 

£ f £ v dv 
[ dv 

Kp = 


kf = : 


kr = 


f B v dv 
f n u F v dv 
f F u dv 
JdB„/dTdv 

J Kv l (dB v /dT) du 


(8a) 

(8b) 

(8c) 

(8d) 


For an optically thick flow, if we assume the radiation 
field is that of a blackbody, then K£ = up and up = «r- 
In contrast, for an optically thin flow, since F v oc £„, it 
follows that Kp = K£ ■ However, the assumption that the 
radiation field is a blackbody is dubious in this regime, 
hence it is unclear whether or not K£ and Kp are even 
related. A full frequency-dependent treatment of radi¬ 
ation would require some transition between optically 
thick and thin regimes. Here, the primary interest is 
investigating the effects of radiation forces that develop 
under optically thick conditions, such that the most rel¬ 
evant single value of the opacity is up = kr. For tem¬ 
peratures < 100 K, the frequency dependence of k,, oc v 2 


leads to an approximate dependence kr oc T 2 (Draine 


2011bI, and some recent simulations investigating radi- 


ation m dusty starburst disks have included this depen - 


dence (|Krumholz & Thompson||2012| |Davis et al.||2014|) . 
For the moderate optical depth as would apply in molec- 
ular clouds of surface density ~ lgcm -2 , the opacity 
would increase by less than a factor ~ 2 from the edge 
to the center of the cloud, and for simplicity we adopt a 
constant k. We allow, however, for a range of values of 
r, as described below. 

We consider four series of runs over a wide range of 
GMC conditions, which we summarize in Table |T| Our 
model parameters are the adopted opacity, k, and the 
initial cloud radius and mass, i? c ioud and M c i ou d, re¬ 
spectively. In addition to the basic input parameters, 
for each model Tabic [I] lists the initial gas surface den¬ 
sity, £ = M cloud /(7r^ c 2 loud ), the initial turbulent ve¬ 
locity dispersion, <j, the initial optical depth through 
the cloud, and the initial free-fall time of the cloud, t 
tff = [7r 2 i? 3 loud /(8GM c ioud)] 1 ^ 2 , We run each simulation 
for 8 times the initial free-fall time, i.e., to tfmai = 8 tg. 

In the K series, we vary the opacity while holding 
-Rcioud = 10 pc and M c i ou( j = 10 6 M 0 fixed, which al¬ 
lows t to var y wh ile £ and a do not. As we shall dis¬ 
cuss in Sect ion pO} simple scaling arguments suggest that 
the most important parameter in determining the fate 
of a cloud with feedback dominated by IR radiation is 
fsdd,* = k4//(47tcG); the K series allows us to focus on 
this dependence. For gas-to-dust ratio ~ 100, realistic 
values of the Rosseland mean opacity for reprocessed IR 
radiati on are likely in the range k ~ 1 — 5 cm 2 g -1 (Se¬ 


menov et al. 2003), but we include a much larger range 
k = 1 — 40 cm 2 g” 1 in the K series to explore the phys¬ 
ical dependence of the results on n. The values at the 
upper end of the opacity range in our models may apply 
for ISM regions that have high dust abundances. Poten¬ 
tially, GMCs with high SFE may have dust abundances 
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Figure 2. Evolutionary history of the fiducial model, with ac = 
20 cm ^_g ~~ 1 , -ftcioud — 10 pc, and M c j ou d = 10 6 M@ (see K20 in 
Table fTF. Shown are the gas mass within the domain (including 
the diffuse background gas, which is initially set to 1% of the initial 
cloud density), the mass that has been ejected from the box, and 
the total mass in star particles. Time is in units of the initial 
free-fall time within the cloud. 

enhanced by self-enrichment. As discussed in Section [l] 
with r > 1 in all runs we expect the reprocessed IR ra¬ 
diation force to dominate over the direct radiation force. 
Also, since er > 16 km s^ 1 , which exceeds the maximum 
expansion velocity of ~ 10 km s^ 1 for H II regions, we 
would not expect significant contribution from H II pres¬ 
sure for clouds in the regime under investigation. 

In the R series, we vary i? c ioud while holding M c i oud = 
10 6 Mg fixed, while in the M series, we vary the cloud 
mass, Mcioud, while holding the cloud radius -ffcloud = 10 
pc fixed. In both the R and M series, we hold k fixed, 
but E, a, and r all depend on both the cloud radius and 
mass, hence vary with either of them, albeit in different 
ways. Finally, in the RM series, we vary both R c ioud and 
Md ou d together such that M c i ou d oc I? 2 loud , while holding 
k fixed. Thus, in the RM series the gravitational poten¬ 
tial well depth and a vary while E and r do not. For 
the R, M, and RM series, we set k = 20 cm 2 g -1 . This 
value is likely larger than realistic values for Solar metal- 
licity. However, we find from the K series that radiation 
does not have significant effects for lower values of k. and 
in the interest of exploring how the dynamical outcomes 
depend on cloud size and mass we must choose a larger 
value of k. 

For computational expediency, we use c = 5 (<j+c s )t in 
runs K30, K40, R7.1, and M2, instead of the larger value 
c = 10(cr + c s )t used in all other runs. Our studies have 
shown that the typical values of T go % are in fact con¬ 
verged with respect to this somewhat relaxed criterion 
(see Figure]]]). 

3. RESULTS 

3.1. Evolution of a Fiducial Model 

We begin by describing the overall evolution of a fidu¬ 
cial case. To illustrate key features of the interaction 
between the turbulent gas of the cloud and the diffuse ra¬ 
diation field that permeates it, we select the model with 
k = 20 cm 2 g _1 , Rdoud = 10 pc, and M c i oud = 10 6 M Q . 
This run has initial properties as listed, e.g., under K20 
in Table IT| Figure [2] shows the history of the mass of gas 
within the box, the mass ejected, and the total mass of 
all star particles. 

Over time, the gas mass steadily declines (slowing after 
2 iff), the ejected mass steadily increases, and the mass 


Figure 3. Effect of resolution on the mass history of the fiducial 
model (see K20 in Table [l|. After 5 free-fall times, the variation of 
the final stellar mass is oforder 5%, indicating that the evolution of 
the simulations is acceptably converged at a resolution of N = 256. 


© 




o 




Figure 4. Effect of the detailed initial turbulent velocity field 
on the mass history of the fiducial model (see K20 in Table [l]l. 
Red, green, and blue curves represent three distinct realizations of 
the initial turbulence spectrum (i.e., produced by three differently 
seeded random number sequences), all with the same initial kinetic 
energy. After 8 free-fall times, the variation of the final stellar and 
ejected masses is of order 10%. 


in star particles increases rapidly between ~ 1 — 2 %, 
and then reaches a plateau. Accretion onto the star par¬ 
ticles effectively ceases after ~ 3 tg, because the strong 
outward radiation force exceeds the inward gravitational 
force. Over the duration of the simulations, ~ 60% of 
the initial mass in the box (or 70% of the initial cloud 
mass) is accreted onto the star particles. Although a 
small amount of gas (~ 10% of the total) is ejected from 
the box at early times (because a portion of the turbulent 
cloud had sufficiently large initial velocities to escape), 
the gas ejection at t x) 3 tg is driven by the radiation 
force. 

The standard resolution of our simulations, N = 256, 
is dictated by the combined constraints of computational 
cost steeply increasing with resolution, and the desire to 
explore a range of parameters. Although the details of 
turbulence models are always subject to resolution, we 
have confirmed that the overall evolution of our simula¬ 
tions is acceptably converged. Figure [3] shows the mass 
histories analogous to Figure[2] comparing results for res¬ 
olutions N = 128 and N = 256 for the first 5 free-fall 
times. Evidently, the variations are only of order 5%. We 
note that different seeds for the random number genera¬ 
tor used to form the initial turbulent velocity field (with 
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Table 1 

Initial Model Parameters 


Model 

K, 

-^cloud 

-^cloud 

S 

a 

/Edd,* 

T a 

% 

c 


(cm 2 g- 1 ) 

(P c ) 

(10 s M 0 ) 

(g cm -2 ) 

(km s -1 ) 



(Myr) 

(10 3 km s _1 ) 

K1 

1 

10 

1 

0.67 

23 

0.068 

1 

0.54 

0.25 

K5 

5 

10 

1 

0.67 

23 

0.34 

5 

0.54 

1.2 

K10 

10 

10 

1 

0.67 

23 

0.68 

10 

0.54 

2.5 

K20 b 

20 

10 

1 

0.67 

23 

1.4 

20 

0.54 

4.9 

K30 

30 

10 

1 

0.67 

23 

2.0 

30 

0.54 

4.9 

K40 

40 

10 

1 

0.67 

23 

2.7 

40 

0.54 

4.9 

R7.1 

20 

10/V2 

1 

1.30 

27 

1.4 

40 

0.32 

5.8 

RIO b 

20 

10 

1 

0.67 

23 

1.4 

20 

0.54 

4.9 

R14.1 

20 

10 V 2 

1 

0.33 

19 

1.4 

10 

0.90 

2.1 

MO.5 

20 

10 

0.5 

0.33 

16 

1.4 

10 

0.76 

1.8 

Ml b 

20 

10 

1 

0.67 

23 

1.4 

20 

0.54 

4.9 

M2 

20 

10 

2 

1.30 

32 

1.4 

40 

0.38 

6.8 

R5M0.25 

20 

5 

0.25 

0.67 

16 

1.4 

20 

0.38 

3.6 

R7.1M0.5 

20 

10/V2 

0.5 

0.67 

19 

1.4 

20 

0.32 

4.2 

R10M1 b 

20 

10 

1 

0.67 

23 

1.4 

20 

0.54 

4.9 

R14.1M2 

20 

10-\/2 

2 

0.67 

27 

1.4 

20 

0.90 

5.8 

R20M4 

20 

20 

4 

0.67 

32 

1.4 

20 

0.76 

6.8 


a The optical depth through the center of the initial uniform cloud, r = 2^.R c i ou dM c i ou d/(| ^cloud) = I kS - 
b The parameters for these runs (in bold face) are identical and refer to our fiducial model. 



Figure 5. Effect of the source size, R*, on the mass history of 
the fiducial model (see K20 in Tablem). After 5 free-fall times, the 
variation of the final stellar and ejected masses is of order 5% for 
runs with R* a factor of 2 larger and smaller than the fiducial size 
R * = 1 pc. 


a given total kinetic energy) can also lead to of order 
10% variation in the mass history of a cloud, as shown 
in Figure [4j Finally, we note that the physical size of 
the source, i?*, used in Equation |3j leads to of order 5% 
variation in the mass histories for values a factor of 2 
above and below the fiducial source size of R* = 1 pc, as 
shown in Figure [5] 

The action of the radiation force on the gas can be 
seen in snapshots of the structure within the cloud at 
different stages of its evolution, as shown in Figure [6j 
In all of the snapshots, the gas density is highly clumpy 
and filamentary due to the turbulence. The radiation 
energy density is highest immediately surrounding star 
particles, and has other local variations due to interac¬ 
tion with the gas. Overall, the radiation energy density 
decreases outward, as does the gas density. The radiation 
flux points primarily radially away from near the center 
of box, where the most luminous star particles are found, 
except in the immediate vicinity of other star particles. 
Inspection of Figure [6] shows that unlike in implemen¬ 


tations of RHD that adopt flux-limited diffusion (FLD) 
methods, the radiation flux vectors (unsealed by magni¬ 
tude) point independently from the gradient in the ra¬ 
diation energy density. Thus, the radiation force is not 
necessarily co-aligned with the negative gradient of the 
radiation energy density, as it is by construction in FLD- 
based implementations. For comparison, Figure [7] shows 
snapshots of £. = / pdz, the gas surface density in the 


sity 

4 


z-direction at the same times shown in Figure 

Over time, the mean density decreases as gas is ex¬ 
pelled from the volume. The gas remains turbulent 
and filamentary throughout the evolution, and expan¬ 
sion does not lead to the kind of simple shell-like struc¬ 
ture that results for non-turbulent clouds (as se en, e.g ., 
for the simple test shown in Fig. 23 and 26 of S013). 
Figure [8] shows the mass-weighted RMS gas velocity, 
Dims = (2E kin /M ga , s ) 1/2 , where Elkin and M gas are re¬ 
spectively the total kinetic energy and mass of the gas 
remaining on the grid, over the evolution of the fiducial 
model. This velocity declines over ~ tg as the initial tur¬ 
bulence dissipates, rises as radiation feedback becomes 
important, and then declines slowly for the remainder of 
the evolution. 

As discussed in Section [I] in addition to limiting the 
star formation efficiency in a given cloud, radiation feed¬ 
back from clusters can also inject momentum into the 
surrounding ISM. Figure [9] shows the total ejected ki¬ 
netic momentum as a function of time for the fiducial 
model, which is ~ 0.47Af c i ou dcr by the end of the run. For 
the purposes of driving turbulence in the ISM and there¬ 


fore r egulating futu re star formation (Ostriker & Shetty 


2011 hereafter OS11), the most important quantity is 
the ratio of total kinetic momentum ejected to the to¬ 
tal mass in stars formed, p*/M*. For the fiducial model, 
p*/M* = 17 km s -1 ; the ratio p*/(M*tr) is equal to 0.78. 
This is c onsistent with the general expectation discussed 


OS 11 that p*/M* from IR radiation feedback will be 


of order the velocity dispersion in the central star cluster 
that forms. 
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Figure 6. Snapshots of slices through the fiducial model, K20, at successive stages in its evolution. Each plot shows the gas number 
density of hydrogen nuclei n // (logarithmic color scale), radiation energy density (contours), direction of radiation flux (vectors), and star 
particles (spheres with logarithmic color scale indicating mass). The slices are through the most massive star particle in each snapshot, 
at (top left) t = 1 tff and z = 2.8 pc, (top right) t = 3 tg and z = 1.3 pc, (bottom left) t = 5 tg and z = —0.78 pc, and (bottom right) 
t = 8 tfj and z = —2.3 pc. Star particles within A z = ±2 pc of the slice are plotted. The 20 energy density contours are logarithmically 
spaced over the data range in each slice. The color scale for the gas density (blue, top) is in units of cm -3 , and the color scale for the star 
particle mass (red, bottom) is in units of Mg. The slice dimensions are 40 X 40 pc 2 . 


In Sections |3.3| and |3.4| we discuss further how evolu¬ 
tion and outcomes of the simulations vary for different 
model parameters. 

3.2. Radiation Structure in a Fiducial Model 

In this section, w e ag ain employ the fiducial model 
presented in Section 3.1 and now use it to explore the 
cloud’s internal radiation structure as well as the compe¬ 
tition between radiation and gravity. We wish to investi¬ 
gate the radiative structure after stars have formed and 
significant feedback is underway, but before the system 
reaches a limiting state, either by accreting the bulk of 
the mass onto star clusters or ejecting it from the simu¬ 
lation domain. Thus, we have chosen to take snapshots 
of the variables at time t = 3 tg. 

3.2.1. Radiation and Gravity in a Spherical System 

Throughout this section, we compare our results to 
those of an isotropic dusty gas sphere surrounding a sin¬ 
gle, centrally-l ocated massive and luminous cluster (cf. 
Appendix A of OS11) (see also Elmegreen 1983; Scoville 


g^_al.| [200li |Krumholz fe Matzner||2009| | Murray et al.| 
2010p . For this system, in steady state the radiation flux 
is F r .* = L*/(47rr 2 ) = 'F M* / (Anr 2 ), and the gravita¬ 
tional field from the cluster is g* = GM* /r 2 . Within any 
radial shell of mass SM, the ratio of the total radial ra¬ 
diation force nF r ,*6M/c to the total radial gravitational 
force g*SM from the cluster is therefore 


/Edd,* = 


eg* 


K'k 

47tcG 


= 0.68 


10 cm 2 g~ 


( 9 ) 


_1_V 

1700 erg s -1 g -1 J 


Including the self-gravity from the gaseous sphere, the 
ratio of the outward force to the inward force on a radial 
shell at r is 


./Edr] ,spii ( r ) 


/Edd,* 

1 + M gas (r)/M* ’ 


( 10 ) 


where M gas (r) is the mass of the gas sphere interior 
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Figure 7. Same as Figure[fi] but for snapshots of £ 2 , the gas surface density integrated in the ^-direction (logarithmic color scale), along 
with a two-dimensional projection of all star particles (spheres with logarithmic color scale indicating particle mass Mi,*). The snapshots 
are at (top left) t = 1 %, (top right) t = 3 %, (bottom left) t = 5 tff, and (bottom right) t = 8 %. The color scale for the gas surface 
density (spectrum, top) is in units of g cm -2 , and the color scale for the star particle mass (violet, bottom) is in units of Mq. The image 
dimensions are 40 X 40 pc 2 . 




Figure 8. Time evolution of the root mean-square (RMS) gas 
velocity in units of the initial turbulent velocity, cr, for the fiducial 
model (K20). 

to r. For any gas mass M &as (r) less than M max = 
A/* (/*Edd,* ^ 1), the outward force on a shell at r ex¬ 
ceeds the inward force. For the fiducial model, with 
k — 20 cm 2 g _1 , /Edd,* = 1-4- As a consequence, if the 
cloud were able to remain spherically symmetric while 


Figure 9. Total ejected radial kinetic momentum in the fidu¬ 
cial model (K20). Momentum is measured in units of pturb,init = 
-^cloud* 7 - 

forming stars at its center, gas out to a radius for which 
A/gas (r)/M* = 0.4 would have the local radiation force 
exceeding gravity. 

If the net efficiency of star formation compared to the 
initial cloud mass is e* = M„/M c i ou d, then the total re- 
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maining gas has M gas (rm ax )/M* = (1 - e*)/e*. Sub¬ 
stitution in Equation ( p~0] ) shows that the local Edding¬ 
ton ratio /Edd,sph(f) would exceed unity at all radii in a 
spherical cloud when a fraction 

(ID 

K \ _1 / 'F \ _1 


£*,sph 


— /Edd,* 


= 1.5 


v 10 cm 2 g 


—i 


\1700 erg s _1 g 


—l 


of the original cloud has been converted to starsjj For 
an idealized spherical system with k = 20 cm 2 g _ , 74% 
of the original cloud would have to be converted to stars 
for the remainder to be forced outward by radiation. Of 
course, gas pressure forces can become important, so that 
gas originally at small radii where /Edd,sph(f) ^ 1 can 
sweep up and expel fluid elements that were originally 
at large radii; this reduces the actual SFE £* i fi na i below 

£*,sph- 

The above discussion suggests that unless k is rel¬ 
atively large, radiation forces will not be able to ex¬ 
pel a substantial portion of a cloud’s mass. In fact, 
larger values of k than would be realistic for IR radi¬ 
ation at Milky Way dust abundances may be required 
for any gas expulsion, even in the limit of highly effi¬ 
cient star formation. For the idealized spherical case, 
the (stellar-plus-gas) gravitational force exceeds the radi¬ 
ation force at the outer edge of the cloud even for e* —► 1 
if /Edd,* < 1- Moreover, for gas very near any star clus¬ 
ter, the dominant forces are radiation and gravity, so 
only if /Edd,* > 1 can continued accretion be halted. 
The condition /Edd,* > 1 translates to k > K cr ;t where 


^crit — 


47tcG 


= 15 cm 2 g -1 


( 12 ) 




1700 erg s -1 g~ 


Although the true structure is non-spherical, this ex¬ 
plains why we have selected the k = 20 cm 2 g _1 model 
for detailed presentation of internal radiation structure: 
for lower-re; models, radiation forces are not expected to 
be strong enough to prevent continuing accretion and 
expel substantial amounts of mass. 

As we shall next discuss, the simple spherical analy¬ 
sis provides useful guidance regarding the relative im¬ 
portance of radiation and gravity, but the solution of 
the full RHD equations for a turbulent medium is more 
complex. In particular, neither the radiation field nor the 
gravitational field is spherically symmetric, and the gas 
density on which these fields act is highly nonuniform. 
The radiation and gravitational forces— locally, averaged 
over spherical shells, or integrated over the whole cloud 
therefore can differ substantially from estimates based on 
simple spherical models. 

3.2.2. Analysis of Radiation and Gravity in the Turbulent 

Cloud 

6 If all the remaining gas is collected in a sing le t hin shell, as for 
the idealized system described in Appendix A of |OSll| self-gravity 
is diluted so that the ratio of radiation to total gravitational forces 
is /Edd = /Edd,*/(1 + 0.5M s h e ii/A^*)• Shell expulsion would com¬ 
mence when a fraction e* idl = (2/p; dd „ — 1) — 1 of the gas is con¬ 
verted to stars; for the case re = 20 cm 2 g —1 this would correspond 
to 56%. 



Figure 10. Angle-averaged radial radiation force per unit mass, 
(re./y/c) (solid), and the angle-averaged total gravitational force 
per unit mass, (chi’) (dash-dotted), both compensated by 4nr 2 , 
for the fiducial run (model K20) at time t = 3 tg. For compari¬ 
son, we plot reL*/c (dashed) and 4irGM t , (dotted), the radiation 
and gravitational forces that would be produced by a point-source 
cluster with the same total stellar mass as is present at t = 3 tg in 
the simulation. The measured radiation force reaches 90% of the 
point-source value at ~10 pc from the center of mass, comparable 
to the radius of the original cloud. At r > 15 pc, 4irr 2 (F r ) asymp¬ 
totes to L*. The measured gravitational force exceeds the stellar 
point-source value at large radii because a significant portion of 
the cloud’s gas has neither been accreted nor expelled at this time. 

To quantify the internal cloud structure in a simple 
way, we compute angle-averaged radial profiles, denoted 
by (•) 4 tt. To do this, we first compute the center of mass 
of the stars, rgM, then linearly interpolate the Cartesian 
grid data onto an N r x x Ng spherical grid with N r = 
Ng = N and N^ = 2N zones centered at tcm, with r £ 
[0, 2i?. c ioud — fCM]Q Finally, we average the interpolated 
data on this grid over all solid angles to obtain a radial 
profile. 

Figure [To| compares, for t = 3 tg, the profiles of the 
two dominant competing effects in the problem: the out¬ 
ward radial radiation force and the inward gravitational 
force. These two forces both scale with the mass (or, 
per unit volume, the density), which we omit here, and 
both fall off as inverse square laws away from the center 
of mass of the system; thus, we compensate by a factor 
of 47rr 2 in each profile. For comparison, we also plot the 
radiation and stellar gravity forces/mass that would ap¬ 
ply in the limit of a single point mass cluster at rcM> 
which appear here as constants (kL*/c and AirGAf*, re¬ 
spectively) in our radially-compensated plot. Far from 
the center of mass of the system, the profile of the radi¬ 
ation force approaches that of a point source, reaching 
90% of this value by r ~ 10 pc, a distance comparable 
to the original cloud radius. However, the gravitational 
force at this time still contains a significant contribution 
from the remaining gas in the system (at this time, the 
total mass in gas is 3.4 x 10 5 M 0 , while the mass in 
stars is 6.6 x 10 5 Mg), so it is slightly higher than the 
point-source cluster limit at large distance. The mean 
specific radiation force exceeds the mean specific gravi¬ 
tational force for most of the profile, indicating that in 
most zones the gas feels a net outward force at this point. 
Figure [2] shows that the accretion rate onto the sink par¬ 
ticles slows between t ~ 2 — 3 tg and essentially stops 

7 Note that rcM -C 2i? c i oud typically, so that 2it c i oud - rcM ~ 

^box/2. 
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Figure 11. Angle-averaged profiles of £ (top) and P rr (bottom) 
for model K20 at time t = 3 tff (solid curves). For comparison, 
we also plot the semi-analytic solutions (dashed) derived from the 
Mi closure for a spherically symmetric system with density and 
radiation flux profiles given by angle-averages of the simulation 
data at this time. The measured data and semi-analytic solutions 
are comparable, showing that the radiation is in a quasi-steady 
state at this time. 


after t ~ 3 tg. We note that in the simple spherical 
model described above, Equation Il2] suggests that radi¬ 
ation forces would exceed gravity forces everywhere only 
when e* —► 0.7 for k = 20 cm 2 g -1 . In fact, e* is only 
^0.57atf~3fff, somewhat lower than this value. 

In Figure B we plot the radiation energy density, 
and the radial component of the radiation pressure 
tensor, (. P rr )ini at time t = 3 tff. For comparison, we also 
plot the semi-analytic solutions for £ and P rr , respec¬ 
tively, derived from the M\ closure for a spherically sym- 
m etric sy stem in steady state (see Equations 110 and 111 


metric sys 

S()l.3j. 


This requires the solution of an ordinary differ¬ 
ential equation (ODE^for the reduced flux, / = F/(c£), 
as a function of the radial profiles (p)^ and (F r ) 4 W . 
The angle-averaged profiles and semi-analytic models are 
comparable for both £ and P rr , indicating that the radi¬ 
ation field is, in an averaged sense, close to quasi-steady 
state for the density distribution at this time. 

Figure 12 shows (again for model K20 at t = 3 tg) the 
radial profile of the Eddington factor, defined by 


/Edd(r) 


(f)Kp r j c) It, 
(pd r $)4n 


(13) 


where the angular averages in Equation (13) are over lo¬ 
cal spherical shells. The quantity /Edd(/]nieasures the 


8 This ODE has a regular singularity at / = 2\/3/5 ~ 0.7 that 
poses numerical difficulties in its solution, resulting in a slight s hift 
in the comparison profile near ~ 12 pc for model K20 in Figure [lT] 



Figure 12. Radial profile of the Eddington fa ctor , f^dd (solid), 
averaged over local spherical shells (see Equation |13| > for model K20 
at t = 3 %. The local Eddington ratio is less than 1 (dotted) over 
most of the domain, indicating that the majority of the gas in the 
cloud feels a net inward force. The ratio of angle-averaged specific 
forces, /Edd,spec = (i^F r / c) / (d r ^) 4 n is also shown (dash-dotted) 
for comparison; this exceeds un ity for most radii. The dashed line 
shows /Edd,* ( see Equation 1 10||, the Eddington ratio for a simple 
spherical model with negligible gas self-gravity. 



Figure 13. Radial profile of xp,F r > the correlation fraction be¬ 
tween p and F r (see Equation |14[ ) for model K20 at t = 3 %. 
Values of Xp,F r less than unity show that p and F r are somewhat 
anti-correlated over most of the domain. 


ratio of the total radial radiation force acting on the gas 
to the total gravitational force acting on that gas, within 
a shell at a given radius. For comparison, we show the 
ratio of angle-averaged specific forces /Edd specMi he., 
with the density omi tted from Equation (13). As previ¬ 
ously seen in Figure [lOj where the angle-averaged spe¬ 
cific radiation force exceeds that of gravity at most radii, 
/Edd,spec( t*) > 1 almost everywhere. In contrast, the ra¬ 
tio of angle-averaged total forces /Edd(0 is below 1 for 
most of the profile. For reference, we also show the quan¬ 
tity /Edd,* defined in Equation (101; this is well above 
both /sdd, spec (r) and / E dd(r) because self-gravity adds 
appreciably to the gravity of the stars at this time. The 
difference between /Edd and /Edd, spec suggests that the 
radiation flux and gas density must be anti-correlated to 
some extent, and highlights the importance of conduct¬ 
ing fully th ree- dimensional simulations. 

In Figure fl3l we plot the profile of Xp,f t - the correla¬ 
tion fraction of p and F r , defined as 


_ (pPr) 4-7T 

xp ' Fr = (pUAFvW 


(14) 


When Xp,F r = 1, the gas density and radiation flux in a 
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Figure 14. Radial profi le o f the cumulative Eddington ratio, 
/Edd,cum( r ) ( see Equation |15| ) , measured inward from the largest 
radius in the interpolation grid for model K20 at t = 3 Since 
/Edd cum( r ) < 1 over the entire domain, the gas feels a net inward 
force. 

given radial shell have uncorrelated fluctuations, whereas 
Xp,F r < 1 implies that th e flu ctuations about the mean 
are anti-correlated. Figure 13 indicates that there is some 
degree of anti-correlation between p and F r , which ex¬ 
plains why /e dd lies below / E dd ,spec 

Figure 14 shows the radial profile of /Edd,cum, the cu¬ 
mulative mfegral of /Edd defined as 


/Edd,cum (l") — 


/ r rmax {pKF r /c) 4 n 47rr 2 dr 
fr ™“ (pd r $)4ir 4-Kr 2 dr 


(15) 


Note that in Equation (151), the volume integrals are cu¬ 
mulative, starting from the maximum radius of the spher¬ 
ical interpolation grid, instead of only over local shells as 
in Equation (13). Once again, this shows that the net 
force that the majority of the gas feels is inward, since 
/Edd,cum < 1 throughout the domain. We note, however, 
that in spite of the net inward force, and the overall anti¬ 
correlation in the fluctuations of density and radiation 
flux, individual fluid elements can experience a radiation 
force (potentially aided by a pressure force) that exceeds 
the gravitational force. This is why, as shown in Fig¬ 
ures [2] and [9l mass and momentum can be ejected by the 
action of radiation forces. 

In analytic models, effects of reprocessed IR radiation 
are sometimes described in terms of a radiation “trap¬ 
ping factor” (e.g., Krumholz & Matzner 2009). This 
measures the enhancements the radiatidnTorce due to 
reprocessed radiation energy trapped in an opaque cloud. 
The trapping factor is defined as a ratio of the radial ra¬ 
diation force on the gas to the force L*/c that would 
apply in the single-absorption limit j£] In Figure 15 


we 


plot the cumulative profile of the trapping factor denned 


by 


f 

/trap,cum (e) = — 


1 ( pnF r / c) 4 ,r 47 tt 2 dr 


L*/t 


(16) 


in comparison to the cumulative optical depth r(r) de¬ 
fined by 

pT max 

r(r) = / (pK) i7I dr. (17) 

J r 

Both quantities are measured radially-inwarcl from the 

9 Note that we do not include UV radiation in these simulations, 
so the trapping factor here only includes effects of IR radiation. 



Figure 15. Radial profi le o f the cumulative trapping factor, 
/trap,cum(r) (see Equation |16| ) , measured inward from the largest 
radius in the interpolation grid (solid), for model K20 at /, = 3 /rr . 
Also shown (dashed) is the optical depth r(r) (see Equation |17| , 
also measured radially-inward. For a simple spherical system with 
a central point source, these would be identical; the discrepancy 
here is explained by the anti-correlation of p and F r and the dis¬ 
tribution of star formation. 

maximum radius of the interpolation grid. For a spheri¬ 
cal, isotropic radiation field of a point source and a spher¬ 
ical gas density profile, /trap cum(f) and r(r) would be 
identical. However, Figure [l}| shows that /trap,cum (f) 
is substantially less than r(rj (and is less than 1 over 
most of the domain). Once again, much of this differ¬ 
ence is explained by the anti-c orre lation between p and 
F r , as demonstrated in Figure [13] In addition, the fact 
that sources are distributed rather than concentrated in 
a single central point implies that the flux is not radially 
directed at small r; the high density in this region there¬ 
fore leads to a greater contribution to r (Equation 17) 
than to /trap,cum (Equation 16). We note that for tire 
present models, the source distribution is in part due 
to the spatial separation of multiple sink particles, and 
in part due to the finite size of each radiation source 
(~ R* = 1 pc; see Equation [3]) . 

3.3. Effects of Varying Opacity 

Here we present the results from the K series, in which 
we vary the opacity while holding the cloud mass and 
radius fixed (see Tab le IT] for input parameter values). As 
explained in Section |2(2| we use a wide range of k; this 
range contains values that substantially exceed the real¬ 
istic mean opacity of dust to IR in star-forming GMCs 
for Solar neighborhood dust abundance, but it allows us 
to study the physical dependence of the cloud outcomes 
on this principal pa ram eter. In particular, the discussion 
leading to Equation [13] suggests that the potential for re¬ 
processed radiation to drive substantial mass loss from a 
cloud depends primarily on k. In high dust abundance 
systems (either galaxies/ISM regions with high metallic- 
ity, or locally dust-enriched individual clouds), n may be 
in the regime where reprocessed radiation forces become 
quite important. 

The comparative properties for Series K are summa¬ 
rized in Table [2j where e* = Af*/Af to t and e gas = 
Mga S /M t ot are measured over the entire computational 
domain, with the total mass defined as Af to t = Af* + 

Afgas + Afej; /Edd,cum(F), /tra P ,cum(r), and r(r) are mea¬ 
sured over the entire spherical interpolation domain for 


r —> 0 (see Equations! 


; spnencal i 

[i5)[EKh|; 


and ( Xp,F r )r is averaged 
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over the spherical interpolation domain. The quantity 
e e j = Mgj/Mtot is measured from the time integral of the 
mass flux through the surfaces of the computational box. 
We report results from analysis at time t = 3 tg, which is 
after significant star formation feedback has begun, but 
before that feedback can destroy the cloud. 

As expected (cf. Equation 12), the star formation effi¬ 
ciency e* decreases as k increases. Correspondingly, both 
the fraction of gas remaining in the computational do¬ 
main, £ gas , and the fraction ejected^ £ e j, increase some¬ 
what with increasing opacity. However, the differences in 
these intermediate efficiencies at t = 3 tg are less than at 
late times, because the dynamical effect from radiation 
requires more time to develop fully. 

In addition to the differences of the efficiencies with n, 
Table|2]shows a number of other interesting effects. First, 
the Eddington ratio integrated over the whole (spherical) 
domain, /Edd,cum (f 0), increases with k . However, 

in all cases it remains less than unity. Moreover, the 
correlation between gas and radiation flux, (Xp,F r )r> de¬ 
creases as k increases for the models with /Edd,* > 1- 
That is, cases with increasing potential for radiation to 
overwhelm gravity (larger /Edd,*) in part compensate for 
this with an increasing anti-correlation of gas and ra¬ 
diation. We further find that the cumulative trapping 
factor, /trap,cum(r —> 0), can be far smaller than the in¬ 
tegrated optical depth r(r —► 0) from the center to the 
edge of the cloud; the ratio is a factor of 4 or 5 for the 
/Edd,* > 1 models. Thus, the estimate tL/c, sometimes 
used in subgrid models (within galaxy formation simula¬ 
tions) to represent the force from reprocessed radiation, 
can subst ant i ally overestimate the true radiation force. 

Figures 16-18 show the time evolution out to t = 8 tg 
of several diagnostic variables. Figure [16| shows the time 
evolution of the mass fractions, £*, £ gas , and £ e j. On the 
one hand, the time evolution of £ gas is similar for all val¬ 
ues of k, with variations of up to 5-10% over all runs in 
the K series. By time t = 8 tg, less that 10% of the orig¬ 
inal mass of the cloud remains, with almost all of the gas 
in each run either accreted onto star particles or expelled 
from the domain. On the other hand, the evolutions of £* 
and £ e j depend strongly on k; as k increases, £* decreases 
and £ e j increases. Figure [16] shows a clear break between 
models in the group with k = 1, 5, 10, compared to the 
group k = 20, 30, 40. These groups have /Edd,* < 1 and 
/Edd,* > 1, respectively, so the break is consis tent with 
general expectations discussed in Section [3. 2. 1| 

There are additional more subtle effects as well. First, 
note that runs K1 and K5 have almost identical mass 
ejection. This is essentially the same as the mass that 
would be ejected in the absence of radiation, due to a 
small fraction (~ 10%) of the initial mass in the cloud 
being unbound. For k = 10, the effects of radiation be¬ 
gins to be seen in reducing £* and increasing £ e j (starting 
at ~ 3 tg). For the large-K group, the final value of £* 
is reached by t = 4 tg in all cases. However, gas ejec¬ 
tion occurs more rapidly in run K40 than in the K20 and 
K10 runs. For large k, the final values of both £* and £ e j 


1(1 It is not sufficient merely to note that the gas has left the 
computational domain through the outer boundary to conclude 
that it has been driven out as a wind, especially considering the 
“diode” outflow boundary condition we employ prevents gas from 
(re)entering through the same boundary. In principle, some of the 
ejected gas might be able to return. 



Figure 16. Time evolution of the mass fractions, e* = M*/Mtot 
(dashed), e gas = M gas /M to t (solid), and £ e j = M e j/M to t (dash- 
dotted), where Mtot = M * + M gas + M e j for the runs in the K 
series. Low-k models KOI, K05, and K10 have /Edd * < 1, while 
high-tt models K20, K30, K40 have /Edd * > 1. There is a clear 
break in £* and £ e j between these groups. Time is in units of the 
initial free-fall time, tff, for each model. 


approach 0.5. 

Figure 17 shows the time evolution of p r ,ej/{M c \ ou gcr), 


the time-mtegrated radial kinetic momentum ejected 
from the grid (in the center of mass frame of the GMC), 
in units of the characteristic initial turbulent momentum 
of the cloud pturb,init = M c \ on g<j. Similar to the results 
for the ejected mass, there is a clear break between the 
low- and high-K models. In all runs, p rie j increases simi¬ 
larly up to ~ 2 tg, due to the initial turbulence expelling 
some of the gas from the grid. After t = 2 tg, the addi¬ 
tional kinetic momentum ejected over the course of the 
simulation is negligible for runs Kl, K5, and K10, where 
/Edd,* < 1) and increases dramatically with k for runs 
K20, K30, and K40, where /Edd,* > 1- 
Finally, Figure [IS] shows the time evolution of a v ir = 
2-Ekin,gas/Agrav,gas: the total virial parameter for the gas 
on the grid. The simulations begin with a v g = 2, i.e., 
with the cloud in the just-bound state. Then, after a 
decrease up to t ~ ts due to the decay of the initial 
turbulence, the virial parameter increases again once star 
formation begins. For runs Kl, K5, and K10, a v ir settles 
into a roughly steady value between 1 and 1.5, i.e., close 
to virial equilibrium. In contrast, for runs K20, K30, and 
K40, the virial parameters rapidly diverge, with larger- 
k models diverging earlier, as the cloud is disrupted by 
radiation from the newly formed stars and the unaccreted 
gas is ejected from the grid. Similar to the results shown 
in Figure [T7J this suggests that a state change occurs in 
systems with /Edd,* > 1 so that the gas quickly becomes 
unbound and is dispersed back into the diffuse ISM. 


3.4. Cloud Evolution Outcomes: Parameter Study 

In this section, we examine the final outcomes for each 
series in our parameter study. Table [3| su mma rizes the re¬ 
sults from each series described in Section [2721 with model 
parameters given in Table [T] From these data, several 
general trends are evident. First, we examine trends in 
the final outcomes of the K series, in which the opacity 
k is varied independently of the other simulation param¬ 
eters. As discussed in Section |3.3[ the star formation 
efficiency £*,/i raa z decreases with increasing opacity, with 
a clear division between models K10 and K20 near the 
critical value of /Edd,* = 1- This division is also clearly 
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Table 2 

Intermediate Outcomes from Varying the Opacity 


Model 

/Edd,* 

e* 

£gas 

£ ej 

/Edd,cum 

(Xp,F r ) r h 

/trap, cum a 

T a 

KOI 

0.068 

0.67 

0.22 

0.11 

0.032 

0.82 

0.074 

0.24 

K05 

0.34 

0.66 

0.23 

0.11 

0.12 

0.77 

0.25 

1.1 

K10 

0.68 

0.63 

0.25 

0.12 

0.17 

0.79 

0.55 

3.0 

K20 

1.4 

0.57 

0.30 

0.13 

0.43 

0.81 

1.1 

5.2 

K30 

2.0 

0.53 

0.33 

0.14 

0.60 

0.70 

1.2 

4.9 

K40 

2.7 

0.50 

0.33 

0.16 

0.99 

0.63 

1.1 

4.1 


Note. — Intermediate outcomes for the K series, in which the opacity is varied. All results 
are given at time t = 3 tff. 


a Measured over the spherical interpolation grid, for r —> 0. 
b Measured over spherical shells, then radially averaged. 



Figure 17. Time evolution of p r ,ej/pturb init> the time-integrated 
radial component of the kinetic momentum ejected from the grid, 
in the center of mass frame of the star particles, in units of the 
initial turbulent momentum Pturb init = -^cioud^ f° r the runs in 
the K series. Only models with /Edd,* > 1 (K20, K30, K40) have 
substantial momentum loss driven by radiation. 



Figure 18. Time evolution of a v i r = 2F'kin/^grav, the total virial 
parameter for the gas on the grid, for the runs in the K series. Note 
that the clouds with k, > /c C rit rapidly become unbound, within a 
few free-fall times, and do so sooner with increasing k. 


seen in the values of p r?e j, the total radial kinetic momen- 
turn ejected from the simulation grid. Table [3] includes 
results for p r ^ e , normalized in three different ways, and 
for all choices the values from k < 10 models are all 
similar, while the values are much larger for k > 20. 

Second, we examine the simulation outcomes in the 
R, M, and RM series as a function of the cloud surface 
density and Mach number. The opacity is k = 20 for 
all of these models, so that /Edd,* = 1-4. On the one 
hand, the star formation efficiencies are all comparable 
throughout the RM series, where the surface density and 


opacity are held constant as the cloud radius and mass 
are varied. On the other hand, as the cloud radius and 
mass are independently varied in the R and M series, re¬ 
spectively, the cloud surface density changes, and there is 
a much stronger variation in the star formation efficien¬ 
cies in these series. This suggests that—all else being 
equal—clouds with higher surface density eject a smaller 
fraction of their initial mass. Interestingly, the trend for 
£*,final with t is the opposite in varying-E models from 
varying-K models. That is, an increase of k for fixed cloud 
mass and radius decreases e*,final) whereas an increase of 
mass or decrease of radius for fixed k increases £*,fi n ai- 
Also, we note that the pairs (MO.5,Ml) and (R5M0.25, 
R10M1) have the same matched values of the initial ve¬ 
locity dispersion (16 and 23 km s -1 , see Table [lj, but 
£*,final varies more strongly for the M series than Tor the 
RM series. That is, surface density appears to be more 
important than the potential well depth for the final star 
formation efficiency of a cloud. 

We further note that although /Edd,* > 1 distinguishes 
models in which e* is reduced by radiation from those in 
which radiation does not affect £*, the values e*,fi na i in 
Table [3] are not consistent with Equation (12). In par¬ 
ticular, the R, M, and RM series all have jEdd,* = 1.4 
and e* jSp h = 0.74, while Table [3] shows a range of values 
for e^finai- Similarly, e* jSp h decreases from 0.74 to 0.37 
from K20 to K40, but the decrease in e* i fi n ai is fraction¬ 
ally smaller. This shows that a simple “Eddington-type” 
spherical model is inadequate for quantitatively predict¬ 
ing the net SFE in a cloud, as controlled by radiation 
feedback. 

In all series, we see that p r , e j / M* , the kinetic momen¬ 
tum ejected per stellar mass formed, is comparable to 
the initial velocity dispersion a of the clouds, to within a 
factor of 2 for all runs with /Edd,* > 1. The values are be¬ 
tween 10 and 40 km s -1 . Interestingly, just as for e^finab 
the ratio p r>e j/M* appears to depend more strongly on 
the cloud’s surface density than its potential well depth; 
all the models in the RM series have essentially the same 
value of p r , e j/M* even though they have varying poten¬ 
tial well depth (i.e., varying a), while p r ^/M* increases 
toward higher E models in the R and M series. 

In principle, it is possible that some fraction of the gas 
that is ejected from the simulation grid might in reality 
ultimately be able to re-collapse. As a proxy for whether 
or not the gas can re-collapse, we compute the escape 
velocity 

/ 2GM c i oud \ 1/2 


^esc — 


\ -^box ) 


(18) 
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Table 3 

Final Outcomes 


Model 

£*,sph 

£*, final 

^ej, final 

Pr,e j 

Pr,e j 

Pr,e j 

Pr,e j / Af* 

Qvir,4 a 





-^cloud cr 

M*cr 

M e j V e sc 

(km s -1 ) 


K1 

15 

0.87 

0.12 

0.16 

0.18 

1.4 

3.9 

1.3 

K5 

3.0 

0.85 

0.12 

0.16 

0.19 

1.4 

4.1 

1.3 

K10 

1.5 

0.81 

0.15 

0.19 

0.23 

1.3 

4.9 

1.1 

K20 

0.75 

0.61 

0.34 

0.47 

0.78 

1.4 

17 

3.3 

K30 

0.50 

0.54 

0.44 

0.73 

1.3 

1.7 

29 

_b 

K40 

0.38 

0.51 

0.49 

0.91 

1.8 

1.9 

38 

_b 

R7.1 

0.75 

0.70 

0.25 

0.60 

0.87 

2.5 

22 

1.9 

R10 

0.75 

0.61 

0.34 

0.47 

0.78 

1.4 

17 

3.3 

R14.1 

0.75 

0.53 

0.39 

0.33 

0.62 

0.86 

11 

5.2 

M0.5 

0.75 

0.50 

0.44 

0.41 

0.82 

0.95 

12 

2.8 

Ml 

0.75 

0.61 

0.34 

0.47 

0.78 

1.4 

17 

3.3 

M2 

0.75 

0.71 

0.23 

0.52 

0.73 

2.3 

22 

1.6 

R5M0.25 

0.75 

0.58 

0.36 

0.60 

1.0 

1.7 

16 

1.5 

R7.1M0.5 

0.75 

0.62 

0.33 

0.54 

0.86 

1.7 

16 

4.0 

R10M1 

0.75 

0.61 

0.34 

0.47 

0.78 

1.4 

17 

3.3 

R14.1M2 

0.75 

0.61 

0.36 

0.41 

0.68 

1.2 

17 

4.8 

R20M4 

0.75 

0.66 

0.27 

0.35 

0.53 

1.3 

16 

1.7 


Note. — All results are given at time t = 8 tff , except for o v i r) 4. 

a The virial parameter is given at time t = 4 tff. 

b In runs K30 and K40, the cloud has become unbound by t = 4 tff and the gas density is at the floor 

value, such that a v i r) 4 is not longer meaningful. 

where i?box = £box/2 = 2i? c i ou( j is the distance from the 
center to the edge of the computational box. In Table [3j 
we compare p rfi j to the product of the total ejected gas 
mass and this escape speed. We see once again that this 
ratio is correlated with £, the initial surface density of 
the cloud, and only for the lowest value of £ considered 
(0.33 g cm -2 in models MO.5 and R14.1) is this ratio 
less than 1. For all other models, the ratio p r ^/M e jV esc 
is greater than 1, indicating that the bulk of the ejected 
gas is likely unbound. 

4. SUMMARY AND DISCUSSION 

In this work, we present the results of a set of numer¬ 
ical RHD simulations of star-forming, turbulent GMCs. 

We focus on the dynamical effects of reprocessed radi¬ 
ation that originates in massive star clusters, in partic¬ 
ular, on the ability of radiation forces to limit further 
gravitational collapse and destroy the GMC by driving 
gas outward. Our models investigate conditions and pro¬ 
cesses similar to those experienced and created by form¬ 
ing SSCs. Our simulations are conducted using the Hype¬ 
rion extension of the Athena code; Hyperion solves the 
two-moment RHD equations using the M\ closure and 
the RSLA (see S013). To represent the interaction of 
IR radiation with the dusty ISM, we adopt the radia¬ 
tive equilibrium condition, such that the absorption and 
emission of radiation by the fluid are assumed to balance 
identically everywhere. Star clusters are treated in an 
idealized manner, with regions of collapsing gas creating 
sink particles that become local radiation sources. We 
adopt spatially-uniform opacity and an isothermal equa¬ 
tion of state for simplicity, but consider models with a 
wide range of opacity, n = 1 — 40 cm 2 g -1 , to test the 
dependence on this important parameter. We also ex¬ 
plore a range of cloud masses and sizes similar to those 
observed in starburst regions, with initial cloud surface 
density £ = 0.3 —1.3 g cm -2 (or 1600 — 6000 M & pc -2 ). 

All models are initiated with a turbulent velocity disper¬ 


sion, <t, such that kinetic and gravitational energy are 
in equipartition. Table IT] summarizes the input model 
parameters. While idealized in several respects, to our 
knowledge this is the first study that has used self- 
consistent, three-dimensional RHD simulations to inves¬ 
tigate the dynamics of IR radiation feedback in cluster¬ 
forming turbulent GMCs. 

In all of our models, evolution to a final st ate occurs 
over several free-fall times (see, e.g., Figure 161. By 
~ 1 tff, gravitational collapse leads to the formation of 
the first star particles, and by ~ 2 — 4 tff, they gain 
most of the mass they will accrete in their lifetimes. The 
evolution of gas mass is similar in all models, with the 
difference that in models with low opacity, essentially all 
the gas is consumed by star formation, while in models 
with high opacity, accretion is halted by radiative feed¬ 
back. For runs with k = 1,5,10 cm 2 g -1 (/Edd,* < 1), 
accretion onto star particles slows after ~ 2 tg; the sur¬ 
rounding gas maintains a quasi-steady value of the virial 
parameter ~ 1 — 1.5 over the next ~ 6 tff (Fig. [l8]), and 
little additional mass is ejected from the simulation box. 
For the runs with k = 20,30,40 cm 2 g -1 (/Edd,* > 1); 
the virial parameter diverges after ~ 2 — 3 tg, and most 
of the remaining gas is ejected from the box. 

The most important parameter in determining £*,finai> 
the net efficiency of star formation over the lifetime 
of a cloud, is /Edd,* = k4'/(47tcG), where T is the 
mean luminosity-to-mass ratio. When /Edd,* < 1, for 
k < 15 cm 2 g -1 (taking T = 1700 erg s -1 g -1 ; see 
Equation [l0| , star formation efficiencies are high (ex¬ 
ceeding 80%); almost all the mass ejected from the box 
(at t < tff) is due to the initial turbulence rather than 
radiation feedback. When /Edd,* > 1, efficiencies are 
lower and decrease as /Edd,* increases, as shown for the 
K series in Tablets] Of secondary importance is the cloud 
surface density £; models with lower £ have lower e* ! fi na i 
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(see the R and M series in Table [3|. [E Models with high 
k and low E have £*,fi na i ~ 0.5. 

The value of /Edd,* is also the discriminant between 
models that have high vs. low values for the ratio of 
ejected momentum to stellar mass formed, p r ^ e] /M* (see 
Table [3]). However, we note that the “high” values 
(jp r ^/M t ~ 10 — 40 km s -1 ) are still quite low com- 


Krumholz & Thompson 


(2012) found, in their 2D RHD 


back sources (cf. Ostriker & Shetty|2011 

; for supernovae, 

this ratio is ~ 3000 km s -i (see, e.g., 

Kim & Ostriker 


20151. For /Edd,* > 1, we find that the ejected momen- 


tum, p r , e j is of order uM*. We also find that p r>e j/M* 
increases with the cloud surface density E. On average, 
gas that is ejected has speed a few times the escape speed 
of the system. 

Our simulations indicate that a large value of k (> 
10 cm 2 g -1 ) would be required for reprocessed radia¬ 
tion to eject a significant proportion of the mass in 
a cloud. This value is perhaps unphysically large for 
IR opacities at Solar neighborhood ab undance even for 
warm dust (see Semenov et al. 2003 for temperature- 
dependent Rosseland mean opacity], with the impli¬ 
cation that reprocessed radiation would not substan¬ 
tially limit star formation in typical Milky Way GMCs. 
However, for galaxies with metallicities higher than So¬ 
lar, or in dust-enriched regions, such large opacities are 
possible. The minimum k for which feedback limits 
star formation in our turbulent simulations is consis¬ 
tent with the prediction of an extremely simple model: 
an isotropic spherical cloud surrounding a central stel¬ 
lar cluster, in which gravity and radiation are the 
only forc es. For this highly idealized situation, Equa¬ 
tion (fl3| ) shows that the critical opacity for which radia¬ 
tion forces can begin to exceed gravity immediately out¬ 
side the cluster (preventing further accretion) is K cr it = 
15 cm 2 g —1 (\R/1700 erg s -1 g -1 ) -1 . In spite of having a 
similar critical transition, however, the simple spherical 
model does not describe the detailed functional depen¬ 
dence of t he n umerical models on k; for example, unlike 
Equation (12), the K series does not show £*,final oc ft' 1 
at k > K cr i71when f Eddi * > 1). 

One difference between the simple spherical model and 
our turbulent simulations is the anti-correlation between 
the gas de nsity , p, and the radial flux, F r (see Fig¬ 


ures 


12 


and 13). Thus, whereas fy,Menhir) oc /Edd,* oc n 


for tne spherical model (Equation [To] ) , the cumulative 
/Edd,cum M ( see Equation [l5| measured in our simula¬ 
tions is not li nea r in n (see Table [2]). In our fiducial 
model, Figure [l2] shows that while /E dd (r) < 1 at most 
radii, /Edd,spec > 1 at most radii; this is a consequence 
of anti-correlation (see Equation fl3] and 14 ). For the 


K series, the anti-correlation between density and radial 
flux increases at higher values of k, and /Edd,cum ( r —► 0) 
remains less than unity for all values of n. Neverthe¬ 
less, even when the cumulative Eddington factor is less 
than unity, some fluid elements become unbound, reduc¬ 
ing £*,final at high enough n. 

11 We note that the relative sensitivity to k and relative in- 
sensitivity to E for the dynamical response of a turbulent, self- 
gravitating cloud to reprocessed radiation is the opposite of the 
dynamical response to direct UV radiation from embedded clusters, 
which strongly depends on S (Raskutti, Ostriker, Sz Sk i nner 2015, 
in pr eparation; see also |Fall et al.||2010| |Thompson Sz Krumholz] 


simulations of turbulent disks using FLD, that the radi¬ 
ation flux is strongly anti-correlated with the matter dis¬ 
tribution. However, the FLD approximation may con¬ 
tribute in part to this result; it is well-known that ra¬ 
diation in this approximation can easily “leak” around 



tions using both FLD and more accurate variable Ed¬ 
dington tensor (VET) methods. They show that the 
anti-correlation is much stronger using the FLD method 
compared to their VET method, since the former does 
not account for the relative insensitivity of the flux to 
strong density contrasts in filamentary structures when 
the optical depths across their widths are of order unity 
or less. Because the M, approximation we have adopted 
does not resolve angular variations in the intensity of 
radiation and can have difficulty capturing the true ra¬ 
diation field in some situations, it will be important to 
check the results we have obtained with other more accu¬ 
rate (but more computationally expensive) RHD meth¬ 
ods such as VET. 



from reprocessed IR radiation. In these treatments, the 
single-scattering UV force, L*/c, is boosted via pho¬ 
ton trapping by a factor tir, the mean optical depth 
to IR through a cloud. In fact, our c omp arisons between 
/trap,cum and t (defined in Equations 16 and[l7| as shown 
in Figure [l5] and Table [ 2 ] indicate that rL*/c may over¬ 
estimate the true reprocessed radiation force by a fac¬ 
tor of ~ 4 — 5. First, the force may be reduced due to 
anti-correlation of the flux and density, as noted above. 
Second, the assumption of a single, centrally-embedded 
cluster may be too naive. A massive cloud may contain 
several clusters in a distributed configuration, especially 
at early stages before these can merge. Radiation forces 
on gas within a cloud with distributed sources is subject 
to cancellation, and only far from the center of mass of 
the distribution would the radiation field approach that 
of a single concentrated source. Meanwhile, the largest 
contribution to the optical depth may be from the dense 
central region of the cloud. 

In addition to accounting for radiation-matter anti¬ 
correlation and distributed radiation sources, it is also 
important for models that apply IR feedback via a sub¬ 
grid model to ensure that the inward gravitational forces 
are consistent with the imposed outward radiation forces. 
This requires that gravity be spatially and temporally 
well-resolved within any clouds where radiation forces are 
applied. If gravity is softened at small scales, then col¬ 
lapse may not occur as rapidly as it realistically should, 
which would give imposed radiation forces an unphysi¬ 
cal advant age. A situation of th is kind might help ex¬ 
plain why Hopkins et al. (2011) concluded that repro¬ 
cessed radiation could play a dominant role in regulat¬ 
ing star formation (for their “HiZ” model), even though 
their smooth particle hydrodynamics (SPH) simulations 
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adopted a value k = 5 cm 2 g -1 that we found leads to a 
negligible reduction in a cloud’s star formation efficiency 

whether for fully turbulent simulations or an idealized 
spherical system. If limited resolution vitiates the direct 
action of small-scale gravity in star-forming clouds, then 
to avoid a gravity/radiation imbalance it would be nec¬ 
essary to incorporate effects of that gravity as part of a 
subgrid model for effects of radiation. More generally, 
the limited resolution of galaxy formation simulations 
precludes direct simulation of many small-scale physical 
processes, but as feedback is crucial to the control of star 
formation, application of subgrid treatments is unavoid¬ 
able. Spatially resolved direct RHD/MHD simulations 
can be used to identify the most important feedback pro¬ 
cesses, and to design and calibrate subgrid treatments 
that capture these key effects. 

In the present work, we have focused on the effects 
of reprocessed IR radiation, with our simulations sug¬ 
gesting that this form of feedback is unlikely to signif¬ 
icantly reduce star formation within GMCs unless the 
dust abundance and opacity are higher than expected 
for Solar metallicity conditions. Alternatively, since the 
opacity n and luminosity-to-mass ratio only appear as 
a product, a top-heavy stellar mass function could in¬ 
crease /Erin.* , p otent ially leading to a reduction in the 
SFE. Turner et al. (2015) present intriguing evidence 
of both dust self-enrichment and a boosted T in Cloud 
D within NGC 5253. Even so, this cloud has estimated 
SFE ~ 0.6, similar to the values of e*,finai we obtain in 
our high /Edd,* and £ ~ 0.33 g cm' 3 models. 

Finally, we note that the direct UV radiation from 
young, hot clusters has an advantage over IR in that op¬ 
tical depths become large even when clouds have much 
lower column densities. Initial study (Raskutti, Ostriker, 
& Skinner 2015, in preparation) of the effects of (non¬ 
ionizing) UV, using RHD models similar to those of this 
paper, suggests that this direct radiation may be effective 
in limiting the SFE within low surface density GMCs. 
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APPENDIX 


FOURIER TRANSFORM POISSON SOLVER WITH OPEN BOUNDARY CONDITIONS 

Equation rf5J) may be rewritten as a discrete convolution over the simulation domain [0, L x ] x [0 ,Ly] x [0 ,L Z \, divided 
into N x N y N7 equal zones. Letting (a,b,c) and ( a',b',c') represent zone-center integer indices, we may define the 
Green function kernel G(x a , y b , z c \ x a ', y b >, z c >) = G(\x a — x a >\, \y b - y b > |, \z c - z c '|) = G(\a- a!\Ax, \b-U\Ay, \c-d\Az) 
as a symmetric function on an extended domain [— L X ,L X \ x [— L y ,L y ] x [— L Z ,L Z ] or equivalently [— N X ,N X — 1] x 
[— N y , N y — 1] x [— N z , N z — 1]. We also extend p(x a , Ub, z c ) over the larger domain, setting the value to zero for a < 0, 
b < 0, or c < 0. Equation © then becomes 

N x -1 Ny-1 N z -\ 

$(x a ,y b ,z c )=G 

a'——N x b'=—N y c'=—N z 

*G(x a ,y b , z c ] x a ’ ,y b ’, z c >) 

xp(x a ',y b ',z c >) AxAyAz. (Al) 


Taking both Gijk and pij b to be 2N X -. 2 N y -, and 2A" 2 -periodic sequences in the indices i, j , and k, respectively, and 
using the discrete analog of the Fourier Convolution Theorem, it follows from Equation (All that 


*&ijk — 


G 


2N X — 1 2N V -1 2N z -l 

E E E 


(2N x )(2N y )(2N z ) ^ _ o n=Q 

X GlmnPlmn 


x exp 


' ( il 

jm 

kn \1 

\2N X 

2 N v 

2N Z J\ 


(A2) 


where Gimn and pi mn are the respective discrete Fourier transforms (DFTs) of the sequences Gijk and pijk- This 
method is computationally efficient since the DFTs can be computed via FFTs. 

In this work, we compute Gimn directly via the DFT of the periodic sequence 

Gijk = -{[(* mod 2N x )Ax] 2 
+[{j mod 2N v )Ay] 2 

+[(fc mod 2N z )Az] 2 } l ^ 2 , (A3) 

where we set C/ooo = 0 to avoid dividing by zero. The DFT is computed once and stored. An alternative to this 
approach is to use the DFT obtained from the finite-difference approximation to the Laplace equation, 


Gimn — 27T 


cos (ttI/N x ) — 1 


Aa; 2 

cos^m/Ny) — 1 

Ay 2 

cos(7 in/N z ) — 1 
Az 2 


(A4) 


(see the Appendix of Gong & Ostriker 2013), where once again we se t Go nn = 0. The approach we adopt here is 


more accurate for long-range forces compared to the result in Equation (A4), although it does not have a closed-form 


analytic expression, and therefore must be stored in an array of size 8 N x N y N z 














